function x  = tridiagblocksolve(A, b, DIM)

% solves block tridiagonal system

[rows,cols]=size(A);

% Determine if number of rows is a multiple of N
if ~mod(DIM, rows)
    error('A must have an integer number of DIM by DIM blocks');
end

% Get number of blocks
num_blocks = rows/DIM;
range = 1:DIM;

% put diagonal and off diagonal block
F = A((DIM+1):end, DIM*0+range);
D = A(:, DIM*1+range);
E = A(1:(end-DIM), DIM*2+range);
L = zeros(size(E));
U = zeros(size(D));

% compute block decompositions
U(range, range) = D(range, range);
for i = 2:num_blocks
    L((i-2)*DIM+range, range) = U((i-2)*DIM+range, range) \ E((i-2)*DIM+range, range);
    U((i-1)*DIM+range, range) = D((i-1)*DIM+range, range) - F((i-2)*DIM+range, range)*L((i-2)*DIM+range, range);
end

% substitute y
y = zeros(size(b));
y(range, 1) = U(range, range) \ b(range);
for i = 2:num_blocks
    y((i-1)*DIM+range) = U((i-1)*DIM+range, range) \ (b((i-1)*DIM+range) - F((i-2)*DIM+range, range)*y((i-2)*DIM+range, 1));
end

% substitute x
x((num_blocks-1)*DIM+range, 1) = y((num_blocks-1)*DIM+range);
for i = (num_blocks-1):-1:1
    x((i-1)*DIM+range) = y((i-1)*DIM+range) - L((i-1)*DIM+range, range)*x((i-0)*DIM+range);
end
